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Hydrodynamic instabilities in shear flows of cohesive 
granular particles 

Kuniyasu Saitoh,*^ Satoshi Takada,^ and Hisao Hayakawa^ 7 

We extend the dynamic van der Waals model introduced by A. Onuki [Phys. Rev. Lett. 94 , 054501 
(2005)] to the description of cohesive granular flows under a plane shear to study their hydrody¬ 
namic instabilities. Numerically solving the dynamic van der Waals model, we observe various 
heterogeneous structures of the density in steady states, where the viscous heating is balanced 
with the energy dissipation caused by inelastic collisions. Based on the linear stability analysis, we 
find that the spatial structures are determined by the mean volume fraction, the applied shear rate, 
and the inelasticity, where the instability is triggered if the system is thermodynamically unstable, 
i.e. the pressure, p, and the volume fraction, 0, satisfy dp/d(j) < 0. 


1 Introduction 

Because flows of granular materials are ubiquitous in nature, a 
better understanding of their properties is crucially important in 
industry and science^. In contrast to usual fluids, inelastic col¬ 
lisions between granular particles significantly influence the dy¬ 
namics of granular flows^-. As a result, the granular flow de¬ 
pends on both the volume fraction and externally applied force 
such as shear rate, where many continuum models have been 
proposed to describe their anomalous rheology such as kinetic 
theory of granular gases under shear— - —, constitutive relations for 
dense granular flows, i.e. the so-called p-I rheology— - —, revised 
non-local models for slow flows of granular materials^ - —, and 
order-parameter descriptions for multiphase (fluid-solid coexis¬ 
tence) flows of granular particles— - —. The interaction between 
grains is also an important factor in flows of granular materials, 
e.g. discontinuous shear thickening of frictional granular parti¬ 
cles— - —, strong shear resistance of dense cohesive granular par¬ 
ticles 36 ! 37 , jammed regime in the flow curve of dense cohesive 
granular materials — 1 2 2 , and instability of freely falling cohesive 
granular streams—. 

Among such a wide range of theoretical and numerical ap¬ 
proaches for granular flows, kinetic theory is one of the most suc¬ 
cessful methods in describing the hydrodynamics of dry granular 
particles— - —. Though the underlying assumption seems to re¬ 
strict its applicability (e.g. the contact duration should be zero), 
kinetic theory gives quantitatively correct predictions of the trans¬ 
port coefficients for rigid granular particles even for moderately 
dense systems^* However, the basic assumption is violated 
once granular particles are aggregated— - —, which is inavoid- 
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able for a collection of cohesive granular particles. Note that 
the cohesive forces can have two different physical origins, i.e. 
van der Waals forces between microscopic powders or capillary 
forces between wet granular particles 6Q i 61 [j. Recently, we stud¬ 
ied the flows of cohesive granular particles under a plane shear 
using molecular dynamics simulations, where we observed vari¬ 
ous spatial patterns caused by heterogeneous aggregates in steady 
states* Since such a heterogeneity cannot be explained by the 
conventional stability analysis of dry granular flowsi£ - — , it is a 
challenging task to explain the hydrodynamic instabilities in co¬ 
hesive granular flows. 


In this paper, we propose an extended dynamic van der Waals 
model originally proposed by A. Onukl-SS to describe hydro- 
dynamic behaviors of a collection of cohesive granular particles. 
Then, we study hydrodynamic instabilities in shear flows of cohe¬ 
sive granular particles with the aid of the dynamic van der Waals 
model. First, we introduce a continuum model of cohesive gran¬ 
ular particles in Sec. [2] where we modify the dynamic van der 
Waals theor ySLSI to include the energy dissipation caused by in¬ 
elastic collisions between granular particles. Then, we numeri¬ 
cally solve the model under a plane shear in Sec. 0 where we 
use the explicit MacCormack schemed for numerical integrations 
and adopt the Lees-Edwards boundary condition— - —. In Sec. [4] 
we analyze the linear stability of homogeneous state to explain 
observed spatial structures in the presence of a shear rate and in¬ 
elasticity. Finally, we discuss and conclude our results in Secs. [5] 
and [6] respectively. In Appendices, we derive linearized hydrody¬ 
namics for the stability analysis (Appendix® and show our per¬ 
turbative calculations of the eigenvalue problem (Appendix [0. 


t In the following, we refer to the van der Waals interactions as the cohesive forces. 
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2 Model 


written as 


In this section, we introduce a continuum model of cohesive gran¬ 
ular materials, where the dynamic van dev Waals theory for multi¬ 
phase fluids—is extended to include the dissipation of energy. 
First, we show hydrodynamic equations of cohesive granular par¬ 
ticles (Sec. io and explain our model of constitutive relations 
(Sec. 12.2D . Second, transport coefficients in the hydrodynamic 
equations are given by the kinetic theory of granular gases, where 
a dissipation rate is also introduced to represent the effect of 
inelastic collisions (Sec. 12.31) . Third, we nondimensionalize the 
hydrodynamic equations and show their homogeneous solution 
(Sec.[2~4l). 


2.1 Hydrodynamic equations 

Let us introduce hydrodynamic fields as the mass density, p = mn, 
velocity field, Ui, and granular temperature 0, T, where m, n, and 
i = x,y,z are the particle mass, the number density, and the coordi¬ 
nate, respectively 0 . Then, the continuity equation, the equation 
of momentum conservation, and the equation of granular tem¬ 
perature in d m -dimension are given by 
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respectively, where we have used the Einstein convention for the 
subscripts (i.j = x.y.f). On the left-hand-sides of Eqs. ©-©, 
the material derivative is introduced as = <9/<9f+with 
the time derivative, d/dt, and gradient, V*. The last term on the 
right-hand-side of Eq. © represents the energy dissipation in the 
bulk caused by inelastic collisions, where we have introduced a 
dissipation rate as £. 


2.2 Constitutive relations 

Next, we discuss the constitutive relations for the stress tensor, 
&ij, and the heat flux, q t . The stress tensor is divided into the 
viscous and reversible parts as 

Gij = fij - ftij , (4) 

where the viscous part is defined as 

%ij = (Vifij + Vjiii ) + 8ij — — fj^ 

(k = x,y,z) with the shear viscosity, f/, and bulk viscosity, f. In 
the dynamic van der Waals theor y 78 ! 79 , the reversible part can be 


% = (p + pi)8ij+MVinVjn , (6) 

where the static pressure is given by the van der Waals equation 
of state, 

nT 2 ™ 
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with the particle volume, vo, and well-depth of the attractive po¬ 
tential for cohesive granular particles, e. In Eq. ©, the diagonal 
part, pi, and higher order gradient, MV/nV jn, with the coupling 
constant, M, represent the increase of energy due to the exis¬ 
tence of interfaces between two different phases. In this paper, 
we adopt the model used in Refs. 78 i 79 for the diagonal part, i.e. 

p\ = — — |Vzi| 2 — MnV 2 n , (8) 

where the coupling constant is assumed to be proportional to the 
temperature as M = 2d 2 v$T with the particle diameter, d, mea¬ 
sured by the range of square-well potential §. It should be noted 
that the coupling term can be derived from a microscopic model 
for thermodynamic interfaces^, but we phenomenologically use 
this expression, because the microscopic derivation for cohesive 
granular particles, so far, does not exist. 

The heat flux is given by 

qi = -kViT - pV t n , (9) 

where the first term on the right-hand-side represents Fourier’s 
law with the thermal conductivity, k. The second term on the 
right-hand-side of Eq. ©, which does not exist in usual fluids, 
is derived from the kinetic theory of granular gases. The physical 
origin of this term can be explained as follows: Inelastic collisions 
in dense regions decrease the kinetic energy of granular particles 
so that the granular temperature tends to be lower than that in 
dilute regions^ - ^! - — . 

2.3 Transport coefficients and the dissipation rate 

Transport coefficients and the dissipation rate of moderately 
dense dry granular particles are well described by the kinetic the¬ 
ory— - —. However, it is still a challenging task to derive those 
for cohesive granular particles, where our attempt to develop a 
kinetic theory of cohesive granular gases is in progress^. In this 
paper, we only study moderately dense systems, where the mean 
volume fraction of granular particles is much lower than 0.5 (but 
is sufficiently dense to be regarded as a finite density system). In 
addition, we assume that the granular particles are nearly elas¬ 
tic and are driven by a small shear rate to keep the low granu¬ 
lar temperature. We have already confirmed that the transport 
coefficients and the dissipation rate of cohesive granular parti¬ 
cles are well approximated by expanding the interaction range 
of a square-well potential with an inelastic repulsive hard-core, at 


$The granular temperature is defined as T = m((v -u) 2 ) /d m n with the velocity of 
granular particle, v, and the local velocity field, u. 

§ The variables with the tilde denote quantities having the physical dimension, while 
those without the tilde, which will be used later, basically denote dimensionless 
quantities. 


^The complete form of the diagonal part is given by p\ — {( nM' — M)/2}|V i n| 2 — 
nMVjn — nT(Vin)Vi(M/T) with M' = dM/dn, where the surface tension is given by 
g = J^ oo M(dn eq /dr) 2 dr with the equilibrium density profile, n eq (r). If the coefficient 
depends only on the temperature, the diagonal part is reduced to the one used in 
this paper^^. 



least, for nearly elastic dilute granular gases^. Therefore, we use 
the transport coefficients and the dissipation rate derived from 
the kinetic theory of inelastic hard-core potentials, where the di¬ 
ameter, d, represents the interaction range of the square-well po¬ 
tential. 

From the kinetic theory of three-dimensional hard-core granu¬ 
lar gases—, the bulk viscosity, shear viscosity, and thermal con¬ 
ductivity are given by 
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respectively, where the explicit forms of dimensionless functions 
of the volume fraction, 0 = vo n, i.e. /^(0), /rj(0), and are 

listed in Table[U Note that their dependences on the temperature, 
Vf, are identical to those in usual hard-core fluids. The transport 
coefficient for the density gradient in the heat flux, Eq. ([9]), is 
given by 

A = 7-3/2 (13) 

64^/7 im 

with the dimensionless function, / M (0), introduced in Table [l] 

The dissipation rate is simply explained by Haff’s law—: The 
decrease of granular temperature by inelastic collisions is pro¬ 
portional to (1 — e 2 )T with the restitution coefficient of granular 
particles, e. The number of collisions per unit time is roughly esti¬ 
mated as xWnVf, where #(0) is the radial distribution function 
at contact. Then, the decrease of temperature per unit time is 
given by dT/dt ~ -(1 -e 2 )nT 3 / 2 = -nTC^ if we assume #(0) ~ 1. 
More precise calculation by kinetic theory— shows the existence 
of an additional term proportional to the velocity gradient such 
that the total dissipation rate is found to be 

C = Ch + (1 -e 2 )f^)ViUi, (14) 

where the dimensionless function, /^(0), is expressed in Table [l] 
The first term on the right-hand-side corresponds to Haff’s law, 
where its explicit form is given by 

os) 

with the dimensionless coefficient, h\(e), as given in Table [lj 


2.4 Nondimensionalization 

We introduce scaling units of the mass, length, energy, and time 
as the particle mass, m, the particle diameter, d, the well-depth 
of the attractive potential for cohesive granular particles, e, and 
a microscopic time scale, t m = d(m/e ) x / 2 , respectively, so that the 
shear rate, y, is scaled as 

s = t m f. (16) 


Table 1 Dimensionless coefficients and dimensionless functions in the 
transport coefficients, where ^(0) is the radial distribution function at 
contact. Here, we have introduced a scaled pressure and derivative of 
the radial distribution function as p* = p/(<j>0) = 1/(1 - 0) - 0/0 and 
X$ = dx/d(j), respectively. 
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Dimensionless hydrodynamic fields are introduced as the volume 
fraction, 0 = vo n, dimensionless velocity field, Ui = (t m /d)Ui, and 
dimensionless granular temperature, 6 = T/e, respectively. Then, 
the hydrodynamic equations 0])-([3]) are nondimensionalized as 
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dm 


OijViUj-Viqi-^W, (19) 


respectively, where we have used the Einstein convention for the 
dimensionless coordinates (/, j = x,y,z) and have introduced the 
dimensionless material derivative as Qf/Sft = d/dt + utVi with 
d/dt = t m d/dt and V; = JV/. In Table [2] we summarize dimen¬ 
sionless forms of the stress tensor, the heat flux, the transport 
coefficients, and the dissipation rate. 

It is readily found that the dimensionless hydrodynamic equa¬ 
tions (H71)-(H9l) have a homogeneous solution, 0 = 0 O , 6 = 0 O > and 
u = uo = (sy,0,0), corresponding to a uniform shear flow, where 
0o, 6q, and u 0 are a homogeneous volume fraction, homogeneous 
temperature, and uniformly sheared velocity field, respectively. 
From Eq. (H9l ), the homogeneous temperature is found to be 


15^(00) 
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Note that a finite value of the homogeneous temperature repre¬ 
sents the balance between the viscous heating and the dissipation 
of energy, where the dimensionless shear rate and inelasticity are 


1-El] I 3 

























Table 2 Dimensionless forms of the stress tensor, aq = (vo/e)^-, the 
heat flux, q t = ( v 0 t m /£d)qi , the transport coefficients, i.e. 

£ = (v 0 r m /mj 2 )|, 77 = (v 0 t m [md 2 )f\, rc=(t m d)k, and p = (; t m /ed 2 )p , and 
the dissipation rate, £ = 7 m £, where the viscous stress, the reversible 
stress, the static pressure, and the diagonal part of the reversible stress 
are nondimensionalized as t*, = (v 0 /e)f/ 7 -, Kq = (v 0 /e)%, p = (vo/e)p, 
and pi = (v 0 /e)pi, respectively. 
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scaled as s 2 ~ 1 - e 2 H 

3 Numerical simulations 

In this section, we numerically solve the dimensionless hydrody¬ 
namic equations (fT71)-([T9l) under a plane shear. We explain our 
numerical setup in Sec. 13.11 and show our numerical results in 
Sec. ESI 

3.1 Setup 

We prepare a periodic LxLxL cubic box with the dimensionless 
system size, L/d = 50, and divide it into N = 125000 (= 50 3 ) small 
cells with the identical volume, d 3 . Next, we randomly distribute 
the volume fraction, dimensionless temperature, and dimension¬ 
less velocity field in each cell around the homogeneous solution, 
i.e. 0 O > and uo = (sy, 0,0), respectively, where the amplitudes 
of fluctuations are less than 10% of the mean values. Then, the 
explicit MacCormack scheme— is used for numerical integrations 
of the dimensionless hydrodynamic equations (H7l)- (H9l) , where 
the dimensionless time increment is fixed to be A t/t m = 0.1. 

To apply a plane shear to the system, we use the Lees-Edwards 
boundary condition which is originally proposed for molecular 
dynamics simulations— and is extended to the finite-element 
method— and the lattice Boltzmann method—. Figure [T| is a 
sketch of our numerical setup, where the centered cube and gray- 
shaded cubes represent the bulk and copies of the bulk ( image- 
cells ), respectively. In this figure, the isosurface in the bulk cor¬ 
responds to the volume fraction, 0j so = 0.35, where the volume 
fractions in the red and blue sides on the isosurface are lower and 
higher than 0 iso , respectively. Then, we move the upper and lower 
image-cells in the opposite directions along the x-axis so that the 
system is sheared by the scaled shear rate, s = t m y. Note that 
our method is different from the remesh procedure— which corre¬ 
sponds to the Sllod algorithm for molecular dynamics simulations, 


|| Note that the homogeneous temperature is an increasing function of time if there is 
no dissipation of energy, where the increase of temperature is equal to the viscous 
heating, i.e. d9o{t)/dt = s 2 r\Q. 


because the external shear is applied only at the boundaries and 
there is no external force in the bulk. 



Fig. 1 (Color online) A sketch of our numerical setup. The centered 
cube represents the bulk, where the isosurface corresponds to the 
volume fraction, 0 iso = 0.35. The volume fraction in the red (blue) side on 
the isosurface is lower (higher) than 0 iso . The gray-shaded cubes are 
copies of the bulk, i.e. image-cells, moving in the opposite directions 
along the x-axis (indicated by the arrows) to apply a plane shear to the 
bulk. 


3.2 Transient dynamics and steady states 

Depending on the mean volume fraction, 0o, dimensionless shear 
rate, 5, and inelasticity, 1 — e 2 , the system exhibits various tran¬ 
sient dynamics and different spatial structures in steady states. 
Figure [2] displays the time evolution of the isosurface, where the 
shear rate is fixed to be s = 3 x 10 -4 . In this figure, the mean vol¬ 
ume fractions and the inelasticities are given by 0 O = (a) 0.80 iso , 
(b) O.90 iso , and (c) 0 iso , and l-e 2 = (a) 3.5 x 10" 7 , (b) 3.0 x 10" 7 , 
and (c) 2 x 10 -7 , respectively. Initially, the isosurface has a ran¬ 
dom structure in space. As time goes on, the density contrast 
starts to grow and the domains merge with each other to make 
a large cluster. If the mean volume fraction is relatively low, the 
cluster is isolated in the bulk so that we observe a spheroidal or a 
droplet like structure in the steady state (Fig.Efa)). On the other 
hand, if the mean volume fraction is relatively high, the cluster is 
elongated along the x-axis by the external shear and we observe 
either a cylindrical structure (Fig.[2]Cb)) or a plate structure (Fig. 
[2](c)) in the steady state. 
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Fig. 2 (Color online) The time evolution of the isosurface for 0 iso = 0.35, where the systems develop from left to right. The volume fraction in the red 
(blue) side on the isosurface is lower (higher) than 0 iso . The mean volume fractions are given by 0 O = (a) O.80 iso , (b) O.90 iS o, and (c) 0i SO , respectively. 
The dimensionless shear rate is fixed to s = 3 x 10 -4 , while the inelasticities are, respectively, given by 1 -e 2 = (a) 3.5 x 10 —7 , (b) 3.0 x 1(T 7 , and (c) 
2x 1(T 7 . 


We then classify spatial structures of the isosurface based on 
the dimensionless wave number, (k x ,k y ,k z ), for the spatial undu¬ 
lation of the isosurface. For example, k x = 0 if the isosurface is 
homogeneous along the v-axis, while k x = k y = 0 if the isosurface 
is homogeneous along both the x- and y-axes, etc. Clearly, the ho¬ 
mogeneous state is characterized by k x = k y = k z = 0. Figure ^dis¬ 
plays typical structures of the isosurface in steady states, where 
we show (a) a droplet (k x = k y = k z ^ 0), (b) a cylinder {k x = 0, 
k y — k z ^ 0), (c) a plate {k x = k z = 0, k y ± 0), (d) a transverse- 
cylinder {k x = k y ^ 0, k z = 0), and (e) a transverse-plate {k x = k y = 0, 
k z ^ 0) structure. Here, we also introduce another case which 
does not belong to any of them as (f) an irregular pattern. 

Next, we map our numerical results onto phase diagrams of 
the dimensionless shear rate, 5, and inelasticity, 1 — e 2 . Figure |4| 
displays the phase diagrams for various mean volume fractions, 
0o, where both the spheroidal and cylindrical structures ( droplet 
and cylinder ) can be observed in relatively low volume fractions 
(Fig- Eta)), while the plate structures ( plate and transverse-plate ) 
appear in higher volume fractions (Figs. Efb) -(d)). In these fig¬ 
ures, the initial homogeneous state is stable if the applied shear is 
large or the inelasticity is small, where the borders between stable 
and unstable regions are well described by the solid lines obtained 
from our linear stability analysis in the next section (Sec.[4]). If the 
system is in the unstable region far from the solid line, i.e. in the 
highly nonlinear regime, the structure in the steady state tends 
to be irregular and strongly depends on the initial condition, e.g. 
the pluses (+) in Figs. 0(a) and (b). Note that we cannot simulate 


systems with highly inelastic situations, i.e. the parameter sets far 
above the solid lines in Fig.[4|Cs 2 <C 1 — e 2 ), because the decrease 
of temperature is too fast to retain numerical stability. 

The system in the steady state is well sheared even though den¬ 
sity contrast is observed in the bulk. Figure \S\ displays the profiles 
of volume fraction, 0(y), and dimensionless velocity field in the 
sheared direction, u x {y), where we have averaged 0(x,yz) and 
u x (x,y,z ) over the x- and z-directions as 


Hy) 

My) 




// 

// 


<j>(x,y,z)dxdz , 


u x (x,y,z)dxdz , 


( 21 ) 

( 22 ) 


respectively. In this figure, the mean volume fraction, dimension¬ 
less shear rate, and inelasticity are given by 0o = O.90i so — 0.31, 
s = 5 x 10 -4 , and 1 — e 2 = 7 x 10 -7 , respectively, such that we 
observe a plate structure of the isosurface in the steady state 
(Fig. Etc)). As shown in Fig. Eta), the initial homogeneous state, 
0(y) = 0o (the open squares), becomes unstable by shear, where 
the density in the steady state (the open circles) is divided into 
a dense region (0(y) ~ 0.5) and dilute regions (0(y) < 0.2) by 
the interfaces around the dimensionless coordinate, y ~ ±10. As 
shown in Fig.EtW, the dimensionless velocity field well develops 
in the steady state (the open circles), though it deviates from the 
initial linear velocity profile, u x (y) = sy (the open squares). 

In our molecular dynamics simulations of cohesive granular 
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Fig. 3 (Color online) Typical structures of the isosurface in steady 
states, where we classify them as (a) a droplet, (b) a cylinder, (c) a 
plate, (d) a transverse-cylinder, (e) a transverse-plate, and (f) an 
irregular pattern, respectively. 
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Fig. 5 (Color online) Profiles of (a) the volume fraction and (b) the 
dimensionless velocity field in the sheared direction plotted against the 
dimensionless coordinate, y = y/d , where we average 0 and u x over the 
x- and z-directions as Eqs. (2T( and (22), respectively. Here, the open 
squares and the open circles are the results in the initial state (t = 0) and 
the steady state (t = 8000), respectively, where the mean volume 
fraction, dimensionless shear rate, and inelasticity are given by 
0o = O.90i so ~ 0.31, j = 5 x 10 —4 , and l-e 2 = 7 x 10 — 7 , respectively. 



Fig. 4 (Color online) Phase diagrams of the spatial structures in the 
steady states plotted against the dimensionless shear rate, 5, and 
inelasticity, l-e 2 . The mean volume fractions are fixed to 0 O = (a) 
O.80 iso , (b) O.90 iso , (c) 0i SO , and (d) 1.10 iso , respectively. The red (blue) 
region represents that the homogeneous state is unstable (stable). Each 
spatial structure is classified as a homogeneous state (©), a droplet (•), 
a cylinder (■), a plate ( a ), a transverse-cylinder (▼), a transverse-plate 
(♦), or an irregular pattern (+). The solid lines are the results of our 
linear stability analysis, Eq. (37). 


particles^, we observed the corresponding spatial structures to 
those displayed in Fig. [3) where the phase diagrams are qualita¬ 
tively identical to those in Fig. [4) Moreover, the profiles of volume 
fraction and dimensionless velocity field are similar to those in 
Fig- El Though the range of scaled shear rate (1CT 4 < s < 10 -3 ) is 
much smaller than that studied in our molecular dynamics simu¬ 
lations (1CT 4 < smd < 1)—, our continuum model well captures 
the dynamics of cohesive granular particles under a plane shear. 

4 Linear stability analysis 

In this section, we analyze the linear stability of the homogeneous 
state to explain the dependence of observed spatial structures on 
the control parameters, i.e. 0o, s, and 1 — e 2 . 

First, we add small fluctuations, 0, 6, and u = (u x ,u y ,u z ), to 
the homogeneous state as 0 = 0o + 0, 0 = Oo + Q, and u = uo + u, 
respectively. Then, we linearize the dimensionless hydrodynamic 
equations (fl7l) - ([T9l) against the small fluctuations. For example, 
the left-hand-side of the continuity equation (fl7l ) is linearized as 

<90 <90 

-^+u-V0~^+syV x 0. (23) 

Here, the second term on the right-hand-side of Eq. (l23t is the 
result of uo • V0 = syV x (j) which explicitly depends on the coor¬ 
dinate, y. With the aid of Eq. J23l) , we linearize Eqs. (H71)-([T9t 
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+ (t)oVj + T7oV 2 )^ + D 0 V,VA, (26) 

= ( 200 V 2 - ^ 0 ) V z 0 -p 0 V z 0 + b 0 V z VA 

+ A) VyV z u y + ^DqV^ + T 70 V 2 ) z2 z , (27) 


00 

0; 


+ ^yV x 0 — 


+ 


(moV 2 + 0)0^0 + (/c 0 V 2 + d> 0 ) 0 

(05V-y — aV x ) % + (05 V x — aV y ) w-y — aV z w z , 


(28) 


where the coefficients, p#, p e , fj 0 , fty, p e , a> 0 , a> 0 , % ic 0 , Mo, 
and 0 , are listed in Tabled 


Table 3 Coefficients in the linearized hydrodynamic equations I24l-l28l, 
where 770 , £o, and /c 0 are 77 , £, and k, in the homogeneous state, 
respectively. The subscripts (0 and 9) represent their derivatives in the 
homogeneous state, i.e. p 0 = dp/d(J), p e = dp/ 00, pe = dp/99, 

Pe = dp/39, 0)0 = 0o d/ 00, and co 0 = dco/99, where 
co = s 2 p-d m (\)9^l2. 


a 0 = 2p e 9o/(d m (\)o), 

Vo = {l-2/d m )po + $o, 

Po = W0 0 , 

ho = uo/0o, 

ko=2Ko/(d m (j)o), 

jlo = 2po/(d m (l)o), 

Pj =P0/0o, 

Pe = Pe/<\> 0 , 

P(j> = Pq/Q 0 , 

Pe = Pd/<\>o, 

C0(j) = 2(00/ ( d m (J)Q ), 

coo = 2(00/(d m 0o). 

0 = 4 WWn0o), 

a = a o + (l-e 2 )0 o /c(0o). 


Second, we introduce the Fourier transforms as 


0 = 

J <k{t)e lkr dk , 

(29) 

0 = 

J 9k{t)e lk l dk , 

(30) 

Uj = 

i J Ujk(t)e lk r dk , 

(31) 


where i and k = (k x ,k y ,k z ) are the imaginary unit and dimension¬ 
less wave number vector (such that the wave number vector is 
given by k = k/0), respectively The Fourier transform of syV x (j) in 


Eq. (l23t is given by (see Appendix IA.4I and Ref.—) 

syVj = - Jsk x ^e ikT dk . (32) 

Thus, the Fourier transforms of the linearized hydrodynamic 
equations (|24l)- (l28t are written as 


(dt ^ dk^ ^ ~ ’ (33) 

where (fa = ( 0 k, 0 k,w x k, u y k,u z k) T is a transverse vector of the 
Fourier coefficients and j£? is a time-independent 5x5 matrix de¬ 
fined as Eq. (1531 ) in Appendix IA.4I 

Third, we introduce a growth rate of the Fourier coefficients as 
0k(0 «= so that the linearized hydrodynamic equation (l33t is 
reduced to an eigenvalue problem, 


(J ¥ + sk x 


00v 


(pk — A<pk 


(34) 


In Appendix [Bj we perturbatively solve the eigenvalue problem 
J34l ) by expanding the eigenvalues, eigenvectors, and matrix into 
the powers of the wave number, k = |k|. In our perturbative cal¬ 
culations, the shear rate and inelasticity are scaled as s ~ 0(k 2 ) 
and 1 — e 2 ~ 0(k 4 ), respectively, so that the homogeneous tem¬ 
perature, 0o ~ s 2 /(l — e 2 ), remains as finite. Then, we find that 
the eigenvalue for the most unstable mode is given by A = A^ 3 ) 
with 


.1(3)^ k 2 

dmfyof 2 


(35) 


(see Eq. (l99l) in Appendix IB,3D . where we have truncated the 
expansion of A^ 3 ) at k 2 and have introduced a coefficient, / = 
\Z a oP(j) T fioPo • Therefore, the eigenvalue is positive if 


dp 

p * = 3j < °- 


(36) 


i.e. the hydrodynamic instability is triggered if the system is ther¬ 
modynamically unstable. Note that the other factor in Eq. (f35l) 
is negative, — 2K$k 2 / d m (j)of 2 < 0- The neutral curve, i.e. = 0, 
is given by the van der Waals equation of state, Eq. $7$, and the 
homogeneous granular temperature, Eq. (|20P . where the dimen¬ 
sionless critical shear rate for the neutral stability is found to be 


^cr — 


2Kd m <t>l(l - {3*1 (e) + 32} (1 - e 2 ) 


15/, (fo) 


(37) 


The solid lines in the phase diagrams (Fig. [4]) are given by Eq. J37l ) 
which well describe the results of numerical simulations. Note 
that there is no fitting parameter in Eq. J37P . 


Our perturbative calculation also agrees with the numerical so¬ 
lution of the eigenvalue problem, Eq. (l34l) . Figure [ 6 ] is a stability 
diagram plotted against the shear rate, 5 , and inelasticity, 1 — e 2 , 
where the solid line is the neutral curve, Eq. (|37l) , and the open 
circles are numerical results of the critical shear rate. Here, the 
LAPACK subroutines^ are used to numerically solve the eigen¬ 
value problem, Eq. (l34l) , where we confirm a good agreement be¬ 
tween our perturbative calculation and the numerical result. As 
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shown in Fig. [7j we also confirm that Eq. J37l ) well describes nu¬ 
merical results with different mean volume fractions, where the 
unstable region increases with the increase of 0o- 

It should be noted that the second term on the left-hand-side 
of the linearized hydrodynamic equation J33l) . i.e. -sk x d(p^/dk y , 
can be eliminated by introducing the time-dependent wave num¬ 
ber vector as k (t) = (k x ,k y — stk x ,k z ), i.e. the Kelvin mode. In 
this case, however, we cannot use an ordinary procedure for the 
linear stability analysis, where the 5x5 matrix becomes time- 
dependent, so that the eigenvalues also depend on time 

and the eigenvectors have to be constructed of Green’s function 
as <pk(0 = / G(k,k / ,r)<Pk'(0)^k / . For the details of this method for 
dry granular flows, see our previous work in Ref.—. 



0 0.1 0.2 0.3 0.4 0.5 0.6 

S xio -3 

Fig. 6 (Color online) A stability diagram plotted against s and 1 -e 2 , 
where the mean volume fraction is fixed to 0 O = 0i so - The unstable (red) 
and stable (blue) regions are divided by the neutral curve, Eq. (37) (the 
solid line), where we obtain a good agreement with the numerical results 
(the open circles). 



0 0.1 0.2 0.3 0.4 0.5 0.6 

S xlO -3 

Fig. 7 (Color online) The neutral curves with different mean volume 
fractions, 0 O , where all the numerical results (the open symbols) are well 
described by Eq. (37) (the lines). Here, increases from O.80 iso to 
1.20i SO as listed in the legend and indicated by the arrow. 


5 Discussion 

In this paper, we have studied hydrodynamic instabilities in a con¬ 
tinuum model of cohesive granular particles under a plane shear. 
The dynamic van der Waals theory for multiphase fluids 78 i 79 has 
been extended to include the energy dissipation caused by in¬ 
elastic collisions, where the transport coefficients and dissipation 
rate derived from the kinetic theory of three-dimensional inelastic 
hard-core potential— were used. 

We have numerically solved the hydrodynamic equations for 
various values of the control parameters, i.e. 0o, s, and 1 — e 2 , 
where the explicit MacCormack scheme^ was adopted for nu¬ 
merical integrations. To apply a plane shear to the system, the 
Lees-Edwards boundary condition—"— was used, where there 
was no bulk shear in contrast to the remesh procedure^. Then, 
we observed heterogeneous structures of the density field in 
steady states, where a spheroidal or cylindrical structure appeared 
if the mean volume fraction is relatively small, while plate struc¬ 
tures appeared in the systems with higher volume fractions. Note 
that such regular structures can be observed in the vicinity of the 
neutral curve, where we observed various irregular patterns in 
highly nonlinear regimes. All the spatial structures observed in 
our molecular dynamics simulations^ have been reproduced by 
the hydrodynamic equations, where the phase diagrams (Fig. ED 
and the profiles (Fig. [5]) are qualitatively similar to those obtained 
in our previous study. 

To explain the dependence of the spatial structures on the con¬ 
trol parameters, we have analyzed the linear stability of the ho¬ 
mogeneous state, where we perturbatively solved the eigenvalue 
problem for the growth rate of small fluctuations. From our linear 
stability analysis, we have found that the hydrodynamic instabil¬ 
ity is triggered if p§ < 0, i.e. if the system is thermodynamically 
unstable. Then, the boundaries between stable and unstable re¬ 
gions in the phase diagrams (Fig. © were well described by the 
neutral curve, Eq. (|37l) , where we also obtained a good agree¬ 
ment between Eq. (f37t and the numerical result of the eigenvalue 
problem, Eq. d34l ) (Figs. [6] and [7|. 

Though the neutral curve, Eq. ([37]), is given by the stability 
criterion, Eq. ([35]) , the eigenvalue, A^ 3 ) ~ k 2 , is isotropic in the 
Fourier space. In other words, the isotropic eigenvalue cannot 
distinguish the observed spatial structures. On the other hand, we 
also find the anisotropic eigenvalue, A^ 4 ) = se x e y — f)ok 2 (see Eq. 
(11 OOP in Appendix IB.3P . where its stability criterion, Eq. Q107P , 
corresponds to the shear-induced instability for usual (dry) gran¬ 
ular shear flows^&. Therefore, the thermodynamic instability, 
p^ < 0, and the shear-induced instability, Eq. 01Q7P , compete with 
each other, where Eq. A1Q7P also depends on the system size, L, 
through the wave numbers. We find that the isotropic eigenvalue 
is always larger than the anisotropic one, i.e. A^ 3 ) > A^ 4 ), because 
our system size, L = 50 d, is too small to observe the shear-induced 
instability for the range of control parameters studied in this pa¬ 
per. In future, further systematic studies of the pattern selection 
for larger systems will be needed as well as the weakly nonlinear 
analysis for the amplitude equation— 6 . 

It should be noted that the temperature increases with time if 
there is no dissipation of energy. Thus, the homogeneous solution 
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is always linearly stable in the absence of inelastic collisions^. In 
our model, however, the mean temperature converges to a finite 
value in the steady state because the viscous heating is canceled 
by the energy dissipation. Therefore, the hydrodynamic instabil¬ 
ity presented in this paper is one of consequences of the dissipa¬ 
tive nature of granular materials. We also stress that the thermo¬ 
dynamic instability, p§ < 0, can be achieved only if the interac¬ 
tion between the particles is attractive. In addition, the stability 
analyses of dry granular shear flows show that the hydrodynamic 
instability is induced only by the layering mode ( k x = 0), while the 
non-layering mode ( k x ^ 0) is always linearly stable^ -71,73-76 . 
Therefore, spatial undulations in the sheared direction (x-axis), 
e.g. droplets (Fig. [3[a)), transverse-cylinders (Fig. [3[d)), and ir¬ 
regular patterns (Fig.[3£f)), do not exist in dry granular systems. 
Thus, our results are also specific to cohesive granular materials. 

Because we studied moderately dense systems with the mean 
volume fractions around 0i so = 0.35, we have used the trans¬ 
port coefficients and dissipation rate derived from the kinetic 
theory of inelastic hard-core potentials—. This assumption may 
be validated if the externally applied shear rate is so small that 
the granular temperature stays in low values, where the macro¬ 
scopic properties of cohesive granular particles can be approxi¬ 
mated by expanding the interaction range of a square-well po¬ 
tential—. However, the microscopic determinations of realistic 
transport coefficients, the dissipation rate, and the coupling con¬ 
stant, M, of cohesive granular materials are important, where our 
attempt to develop a kinetic theory of cohesive granular gases is 
in progress^. 

Moreover, the effects of gravity and microscopic frictions be¬ 
tween the particles should be examined for practical applications, 
and the influence of the boundary condition, e.g. the study with 
remesh procedure^ or physical boundary conditions, is also im¬ 
portant. 

6 Conclusion 

In conclusion, the extended dynamic van der Waals model can 
describe cohesive granular flows under a plane shear, where the 
hydrodynamic instabilities are well characterized by the neutral 
curve obtained from the linearized hydrodynamics. The various 
spatial structures observed in simulations appear in the unstable 
region, where the hydrodynamic instabilities are triggered if the 
system is thermodynamically unstable, i.e. p^ < 0. 
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A Linearized hydrodynamics 

In this Appendix, we linearize the dimensionless hydrodynamic 
equations (H7l)-(H9l) around the homogeneous solution, 0 = 0 O , 
6 = Qo, and u = uo = (sy,0,0), where 0o, and uo are the homo¬ 
geneous volume fraction, the homogeneous temperature, and the 
linear velocity field, respectively. Here, we add small fluctuations, 
0 , 6 , and u = (%, u y ,u z ), to the homogeneous fields as 0 = 0o + 0, 
6 = 0q + 0, and u = uq + u, respectively. In the following, we lin¬ 
earize the pressure, shear viscosity, and dissipation rate as 


p{<j>,e) 

- Po + p^tjy + ped , 

(38) 

?7(<M) 

^ rio + Vffi + VeQ, 

(39) 

C(M) 

— Co + 0 + Ce + 

(40) 


respectively, where po, Po> and Co are the homogeneous values of 
p, p, and C> respectively, and the derivatives in the homogeneous 
state, i.e. p#, pq, p^, Po, C</>> an d Ce> are listed in Tabled 


Table 4 Dimensionless coefficients in Eqs. (38)-(40). 

P</> = 0 o/ (1 — 0o) 2 - 20o , 

Pd = 0o/(l-0o), 

P<j> = (5v%/160f)^/r](0o)/^0o , 

Pe = (5/32v5f^)/ T? (0o) , 

c 0 = {(3/11 + 32) /24} y/m (1 - e 2 ) (z-+ 00*0) , 

Ce = {(3 h + 32)/48} v / ^b(l - e 2 )0o* , 
cOfj, =s 2 riQ-(d m /2) (Co T 0oC<^>) #o , 

= s 2 Pe ~ (d m /2) (Co + #oCe) 0o , 

X<$> = -367T(7T0 - 15)/(tt0 - 6) 4 , 

v 0 = (7r/5)(l + e)(x + 0X0), 

d fr \/00 = ( h 3 -h 2 )~ l (1 + 2V/3 )dx~ l /d(j) 

+(2/3)4(0)v 0 + (3/45)(32 - hi) (v + 0v 0 ) , 
dx~ l /d(f) - 7T(7T0> - 15)(7T0 -6) 2 / {9(7T0 - 12) 2 } . 


A.l Material derivatives and the continuity equation 

The material derivatives in Eqs. (fl7l) -([T9l) are linearized as 




dX 

~dt 


+ syV x X 


(X = 0, 0 , Uy . U Z ) , 


(41) 


@u x 

Q)t 


dux 

dt 


+ syX x u x + su y , 


(42) 


where the last term on the right-hand-side of Eq. (l42l ) is the result 
of u-V(sy) = su y . Then, the left-hand-sides of Eqs. tT7l -tl9\) are 
linearized as 


^0 





dt 


-+syV x <i> , 

(43) 

(~yf+ syV x u t + SaSUy'j , 

(44) 


(45) 


respectively. 

Because the homogeneous solution is incompressive, i.e. V • 
uq = 0 , the velocity gradient is linearized as V • u ~ V • u so that 
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the linearized continuity equation is given by Eq. (|24l ) . 

A. 2 The equation of motion 

Next, we linearize the dimensionless equation of motion, Eq. 
(fl 8 l) . Since the diagonal part of dimensionless reversible stress 
tensor is linearized as p\ ~ — 20 o 0 o^ 2 0 , the dimensionless re¬ 
versible stress tensor is reduced to Kij ~ poSij + ftij with the first 
order term, 

ftij = {(/ty -200%V 2 ^) (ji+poG^dij , (46) 

where its off-diagonal part is zero. The diagonal part of dimen¬ 
sionless viscous stress tensor is linearized as in ~ in with 

Zii = 2t\qV i&i + - 2-Vo^J V k u k , (47) 

where the dimensionless bulk viscosity, is defined in the ho¬ 
mogeneous state (note that the first term on the right-hand-side, 
VjM;, should not be summed over the subscript, i). The dimen¬ 
sionless shear stress and the other off-diagonal parts of the di¬ 
mensionless viscous stress are linearized as T*y ~ sijo +and 
lij ~ lij (ij ^ xy,yv), respectively, where the first order terms are 
given by 

% ty — t?o (V x w-y + V-yW x )+5 , (r/00+ ?70 0) , (48) 

% = Vo(ViUj + VjUi) (ijj^xy,yx), (49) 

respectively. Then, the dimensionless stress gradient is linearized 
as 

^ 7 j G xj — ^ s V(j>^y ~ (^P(j) ~ 20 o$o^ ^ Vx 0 

+ { sr lo^y ~Pd^x) 0 + + t/oV 2 ^ u x 

■E^oV^V-yM-y + U()V Z V X W Z , 

where we have introduced the sum of dimensionless viscosities as 
Do = (1 —2/d m ) r]o + £o- Therefore, the linearized dimensionless 
equations of motion are found to be Eqs. (l25l) -([27l). 

A. 3 The equation of temperature 

Finally, we linearize the equation of dimensionless temperature, 
Eq. ( fl9l ). It is readily found that On^m ~ —po^k^k> G ij^i u j — 0 
(j j and ij 7 ^ xy,yx), o^V^w-y — sTioV x Uy, and Oyx^yUx — s 2 p o + 
s 2 {r \00 + rjeG) +2srjoV y u x + srjoV x Uy Thus, the dimensionless to¬ 
tal power is linearized as 

G ijVi u j — S 2 P0 + 5 ' 2 ( r /00 + P6@) +2sr/o (VyUx + VxUy) — PO^k^k • 

The heat current and correction term for the dissipation rate are 
linearized as V*# ~ /c o V 2 0 + PoV 2 (j) and -(1 - c 2 )/^(0)V^^ ~ 
— (1 — c 2 )/ ( r( 0 o)V/ : %, respectively, where Kq and qo are the dimen¬ 
sionless thermal conductivity and dimensionless transport coeffi¬ 
cient proportional to the density gradient in the homogeneous 
state, respectively. Thus, the linearized dimensionless equation of 
temperature is given by Eq. (f28l) . 


Note that the zero-th order equation, co 0 = s 2 ri ( j ) - 
(d m /2)0o0oCo = 0, represents the balance between the vis¬ 
cous heating and energy dissipation in the bulk, where the 
dimensionless homogeneous temperature, Go , is given by Eq. 
( 1201 ) . 


A. 4 The Fourier transforms 


We introduce the Fourier transforms of the small fluctuations as 
Eqs. (|29l)-(l3Tl). In the linearized continuity equation (|24l ) , the 
second term on the left-hand-side explicitly depends on the y- 
coordinate, which is transformed as 

syV x $ = ~ Jsk x ^je ikr dk , (50) 

where we have used ^ {d^e lkr /dky)dk y = 0. Therefore, the 
linearized continuity equation in the Fourier space is given by— 

^-sk x -2-^<j> k = <j) 0 k-u k . (51) 

Similarly, we transform the linearized equations (I25l) - (l28l) into 
the Fourier space to find 

{d~t~ skx dk;) (pk= ^ (pk ' (52) 

where we have introduced a vector of the Fourier coefficients as 
<Pk = (0 ki Qk, u xk: u yk: u zk) T and each component of a 5 x 5 matrix, 

(53) 

is defined as 

^\\ = ~^12 = 0 , ~^13 = 0(A , =^14 = 0o ky 5 =^15 = 0O^z j 

J^2i = — fiok 2 , J^22 = (Oq — Kok 2 , J&23 = a k x — bsk y , 

J ^24 = aky — bsk x , J £?25 = , 

-S?31 = sfitpky ~ (p$ + 26 0 k 2 ) fe* , ,S?32 = sfieky - p e k x , 

-£33 = 7 -S?34 = -v 0 k x ky -s , ^35 = -v 0 k z k x , 

-S?41 = sf]qk x - + 200^3) fey , ^?42 = sfjek x -pek y , 

Jz ?43 = -Vok x k y , ^f 44 = -u 0 fey - 7]ofe 2 , -£45 = -v 0 k y k z , 

^51 = - +20 o fe 2 ) k z , i ?52 = -pgk z , ££33 = -v 0 k z k x , 

-S?54 = -Uofeyfez , ^55 = ~%k 2 z - Tfofe 2 , 

with the dimensionless wave number, k= |k|. 
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A. 5 Transverse and longitudinal modes 


We represent the vector, (p^, in a linear combination of the unit 
vectors, 


n l 

= (1,0,0,0,0) T , 

(54) 

n 2 

= (0,1,0,0,0) T , 

(55) 

n 3 

= {0,0,e x ,e y ,e z ) T , 

(56) 

n 4 

= (0,0,-e^e y ,e ± ,-e y e^) T , 

(57) 

n 5 

= (0,0, —e^ ,0,e x ) T , 

(58) 


as (jpk = J ^ =1 ajikj, where the scaled wave numbers are introduced 
as ej = kj/k, ej = kj/k±_ (j = x,y,z), and e±_ = k^/k with k± = (k 2 — 
ky ) 1 / 2 . Here, n 3 is parallel to the dimensionless wave number 
vector, while n 4 and n 5 are perpendicular to it, i.e. n 3 k = k and 
n4 k = n5 • k = 0. 


sk x U~ 

1 (dU / dky) =(P a p), iS 

given by 



P\\ 

H 

K> 

II 

■ 4 ^ 

II 

II 

0 , ^13 = 0o k , 



P21 

= C 0 (j) - iM)k 2 , P22 = 

= a>e ~ kok 2 , 



P23 

= ak — 2 bse x k y , P 2 4 = bsk x (e y ej — ej_) 

, P 25 = 

bsk y e£ , 

P31 

= Isri^exky - (pq, +20o& 2 ) k, P^ = 2 sr)oe x k y - 

~Pek , 

P33 

= - (Oo + Po)k 2 — se x e y , P34 = — 2 se x e± , P35 

= 0 , 

P41 

= sfj 4 >k x (e x ~e y ej-) 

, P 42 = sfl e k x (e ± - 

■e y e y ) , 


P43 

= se x , PiA = se x e y 

- m k 2 , ^45 = 0 , 



P51 

= -sf\$k y e£ , P 52 = 

= - si) g k y e£ , 



P53 

= se y e z ■ p 54 = se z 

> P55 = —fjok 2 ■ 


(63) 


B Perturbation theory 

In this Appendix, we perturbatively solve the eigenvalue problem, 

P+Skx Jk^j ^ = ’ ( 64 ) 

where we have introduced the growth rate which is equivalent to 
the eigenvalue, A, as (fe(t) e^. Here, we also define the left- 
eigenvector, t/ 4 , as 


A new vector, (p k = {a\, is defined as (jpk = U(p^ 
with a matrix, U = ( 111 , 112 , 113 , 114 , 115 ). Then, the linearized hydro- 
dynamic equation (f52l) is transformed asH 

<5,) 

where the left-hand-side is reduced to 


(p T + sk xjj^j$ k=Aftk- (65) 

At first, we introduce a small parameter, e, for the perturbative 
calculations, where the dimensionless wave numbers are scaled 
as 

k= £q , kj_ = £q j_ , k x = £q x , k y = £q y , k z = £q z , (66) 


(l.h.s) 


<9«Pk 

dt 


-ski 


d(pk 

dk v 



du 

dky 


0 k 


with a matrix, 


sk x U~ l 


dU 

dky 


/0 0 0 0 0 \ 

0 0 0 0 0 

0 0 0 —se x e± 0 

0 0 se x e± 0 0 

\0 0 0 0 0 / 


(60) 


(61) 


Therefore, the linearized hydrodynamic equations are rewritten 
as 

(62) 

where each component of the matrix, P = U~ l J?U + 


which do not change the scaled wave numbers, i.e. ej = kj/k = 
qj/q, ej =kj/k ± =qj/q±_ (j=x,y,z), and e± = k ± /k = q ± /q. The 
dimensionless shear rate and inelasticity are scaled as 

s = £ 2 s , l-e 2 = £ 4 g, (67) 

respectively, so that the homogeneous temperature, 6o ~ s 2 /(l — 
e 2 ) = s 2 /q , remains as finite. 

Then, we expand the eigenvalue, right- and left-eigenvectors 
into the powers of £ as 


A 

— + £ 2 A 2 + • • • , 

(68) 

0k 

= <po + £<pi +£ 2 ^2 + ... , 

(69) 


= vb + £Vi+e 2 V2 + --- , 

(70) 


respectively. The matrix, P, is also expanded into the powers of £ 

** The scaled wave numbers satisfy the relations, e 2 x + e 2 + e\ = 1, e 2 + e\ = e]_, e x e x + as 
e 7 e}r = e±, and e x e\ = e,-. Because the unit vectors, n, ( / = 1 ,... , 5 ), are orthonormal, 

the inverse matrix of V is equal to the transposed one, i.e. U~' = U T . P = EqA + t{t /B + sC} + £ i {q i F + i^rG} + ... (71) 
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with the matrices, 


and 


( 0 

0 

00 

0 

o\ 



0 

0 


0 

0 



-Pt 


0 

0 

0 

5 


0 

0 

0 

0 

0 



V o 

0 

0 

0 

(V 



/O 0 

0 



0 

o \ 

0 -i£o 

0 



0 

0 

0 0 

-(Oo + ijo) 


0 

0 

0 0 

0 



-no 

0 

\0 0 

0 



0 

-no) 

/0 0 

0 


0 


0\ 


0 0 

0 


0 


0 


0 0 

~ e x e y 


-1e x e h 

0 


0 0 

e 1 


e x e y 


0 


\0 0 



e z 


<v 



B.l The 1st order equation 



= -j=j(P4,Pe,+if, 0.0), 

(80) 

V / o 2) 

= 0,0), 

(81) 


= -^(ao,-0o,0,0,0) , 

(82) 


= (0,0,0,1,0), 

(83) 


= (0,0,0,0,1), 

(84) 


respectively. 

B.2 The 2nd order equation 

Because the three eigenvalues, A^ (/ = 3,4,5), are degenerated 
to zero, we replace the right-eigenvectors, (p^\ with a linear se¬ 
ries ■■ 

( 85 ) 

where the coefficients, (m = 3,4,5), will be determined in the 

following. Then, the second order equation is found to be 

qA(p P + + sq x -f-Aj (p^ = ^hm <Po m) , (86) 


The first order equation is found to be 


where we have introduced = q 2 B + sC. Multiplying both sides 
by the left-eigenvectors, (n = 3,4,5), we find that the first 

(n) 

term on the left-hand-side vanishes (because of i//g ’A = 0) and 
the second term on the left-hand-side is reduced to 


qAfy o = Ai <po , qA J \\rl = $ , (72) 

where the five eigenvalues are given by 

= xf> = ifq, A, (3) = A, (4) = Aj (5) = 0 (73) 

(note that the superscripts represent different eigenmodes). 
Here, we have introduced a constant as 

f=Jaope + * \ pp + ~T^T • (74) 

v y d m pQ 

The corresponding right- and left-eigenvectors are found to be 


•Po 0 = 

-^-(0O,ao,-i/,O,O) T , 

(75) 

<Po 2) = 

1 x 
-^-(0O,a o ,+//,O,O) , 

(76) 

II 

m 

j{Pe,~P<t>, 0,0,0) T , 

(77) 

II 

(0,0,0,1,0) T , 

(78) 

II 

in' 

(0,0,0,0,1) T , 

(79) 


+ hm<p^ 

(n) yy (m) 7 (/) _ (n) (m) dhrn 

— i//q <Pq + sq x y { o <Po 

dq y 

= <p^ + sq x ^—, (87) 

because (p^ 77 ^ is independent of the wave number and we have 
used * j/q^cPq^ = S nm . Thus, Eq. (l86l) becomes a linear simultane¬ 
ous differential equation for the coefficients, h„ \ i.e. 

_ dh$ (n) yy (m) j (/) ,(/);(/) 

sq x —~. h 9o — ^2 h n , 

dq y 

where the matrix is explicitly given by 

f-gq 2 o o \ 

= 0 se x e y -fi 0 q 2 0 

\ 0 se z -rjoq 2 ) 

with a constant, 

_ 2 / c 0 pq 
8 ~ dmQof 2 


( 88 ) 


(89) 


(90) 
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The differential equation (l88l) for n = 3 is given by 


and 


sqx-^—h^ = (x^ +gq 2 ) , (91) 

where the eigenvalue and solutions for l = 3 are readily found to 

x£ ) = -gq 2 , hf ] = 1, hf ] =hf ] = 0, (92) 

respectively which also satisfy the differential equations (f88t for 
n = 4,5. The other eigenvalues, and xff\ are different from 
X so that = 0 to satisfy the first differential equation 

(f9T1 ) . Then, the differential equations (l88l) for n = 4,5 are given 
by 

* qx l!q- h ^ = ~^ exe y + f l° ? 2 ) /i 4 ) > (93) 

^ qx Jq h ^ = (^2 Z) + W 2 ) ^ -s e z$ 1 (94) 


respectively (/ = 4,5). Here, we choose the eigenvalue for / = 4 as 
A 2 (4) = - ho# 2 t0 find the solution as = 1. Then, the last 

differential equation (1941) is reduced to 



( 4 ) _ ^,( 4 ) 
5 5 


Wx ’ 


(95) 


where its solution is found to be 


J 4 ) = 


tan 



(96) 


Similarly, if we choose the eigenvalue, X^ = ~Poq 2 , the differen¬ 
tial equations (l93t and (f94t are reduced to 



( 5 ) _ 

4 — 


_^J5) 

2 n 4 

q z 


4-h? - --*-»?>, 


dq 


Wx 


(97) 


respectively, where their solutions are readily found to be /z^ = 0 
and /z^ = 1. 


B.3 Unstable mode 


In summary, the eigenvalues and corresponding right- 
eigenvectors are given by 


Ad) = 

-A (2 ) = —ifk , 

(98) 

a® = 

-gk 2 , 

(99) 

A< 4 ) = 

se x e y - ijok 2 , 

(100) 

Ad) = 

1 

-=3i 

0 

(101) 


<Po 1) = 


( 2 ) 

% = 


^ = 


% 


( 4 ) 


-^-(0O,a o ,-;/,O,O) T , 

1 T 

(00 ■ flo, +if, 0,0) , 

y(Pe,-^,0,0,0) T , 


k xk ; tan (-*:)*> 


( 102 ) 

(103) 

(104) 


= [0,0,0,1,^-tan-q-^. 


^ = 


’ V 

(0,0,0,0,1) T , 


(105) 

(106) 


respectively, where these results are consistent with those in 
RefJS. 

Because the first two eigenvalues, A^ 1 ) and X^ 2 \ represent prop¬ 
agating modes and the fifth eigenvalue is negative, A^ 5 ) < 0, the 
third and fourth eigenvalues, Aand X^ 4 \ can be unstable, 
where A^ 3 ) and A^ are isotropic and anisotropic in the Fourier 
space, respectively The isotropic eigenvalue, A^ 3 ), is positive if 
g < 0, i.e. the system is thermodynamically unstable, p§ < 0, which 
could happen in the van der Waals model. On the other hand, the 
anisotropic eigenvalue, X^ 4 \ is positive if se x e y > pok 2 , i.e. 


kxky . f 5v/l5/q(fr)) 3 / 2 ) 1 

k 4 \ l6n^ y /d m (3hi+32)x(M f VT ' 

Eq. J1Q7D corresponds to the shear-induced instability for usual 
(dry) granular shear flowsZ&, where we have used ij 0 = po/(/)o 
and Eq. (l20l) . In cohesive granular shear flows, the thermody¬ 
namic instability, p§ <0, and shear-induced instability, Eq. (U07D , 
compete with each other, where the second criterion, Eq. dl07D , 
also depends on the system size, L, through the wave numbers, 
k x ,k y ,k ~ 1/L. Such a system size dependence of the shear- 
induced instability is consistent with the previous stability analy¬ 
ses of dry granular shear flows— ^2 . Note that the isotropic eigen¬ 
value is always larger than the anisotropic one, i.e. A( 3 ) > X^ 4 \ in 
our system size, L = 50 d, with the range of control parameters 
used in this study. 
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